function [lh,sp]=gensMSM(y,x,sig,b,P,ind,lags)
% Code implements Hamilton's Filtering algorithm
% If ind=1 
% The model is 
% y(t)-M(S(t)) = x(t)b(S(t))+eps(t)sig(S(t)), t=1,...,T
%      x(t) is kx1 
%      ns number of states 
%      b is k x ns
%      sig argument is 1 x ns
%      P(S(t)=j|S(t-1)=i)=pij 
%      The i-th column of P corresponds to P=[pi1 pins]'
%      as in Hamilton (94)  
%      eps(t) i.i.d. N(0,1), independent of current and past x, past y
% 
%  [P(S(0)=i)]=pbar, the steady-state distribution implied by H.
% If x contains lagged y's, 
%   those occurring in x(1) are conditioned on, i.e.
%   treated as non-random.
% Output 
% lh  This function evaluates the likelihood of the sample 
%     (actually, a vector whose product is the likelhihood) 
%      and, for each date t in the sample,
% sp   p(S_t|data_t) as a function of S_t.  
%
ns=size(H,1);
T=size(y,1);
m=size(b,1);
if size(b,2)~=ns, error('b/H dimension mismatch'),end
if size(sig,2)~=ns, error('sig/H dimension mismatch'),end
if size(x,1)~=T, error('x/y length mismatch'),end

% Columns of P must add to one 
[v,d]=eig(P);
% Check that 1 is an eigenvalue of the H matrix i.e. MC
j1=find(abs(1-diag(d))<eps*100);
if isempty(j1) | any( sig<eps*100 )
   lh=-1e80*ones(T,1);sp=.5*ones(T,2);
else
    % Check the ergodicity of the Markov Chain 
    % Find the ergodic probabilities 
    pbar=v(:,j1)/sum(v(:,j1));
    
    if size(pbar,2)>1
        pbar = mean(pbar,2);% in case of non-ergodic H matrix (more than one unit root).
    end
    pst=pbar'; % Row vector
    lh=zeros(T,1); % time series of likelihood increments, like standardized residuals
    sp= zeros(T,ns); % Filtered state probabilities.
    sr2pfac=1/sqrt(2*pi);

    if ind == 2 
        
    
    for it=1:T
      % Hamilton's Algorithm  
      
      %Step 1: Forecast p(S(t)=j |data_t)= j-th column of H'*pst  i.e. 
      %              sum[ ( S(t)=j|S(t-1)=i )* p(S(t-1)=i data_t-1 ) ]
      
      %--------------------------------------------
      % Step 2: Find the joint density f( y(t) , s(t) | y(t-1) )
      % ----------------------------------------------------------
      % sp is then the various components of f( y(t) , s(t) | y(t-1) ) = 
      %     f( y(t)| s(t) ) * p( s(t) |  y(t-1) )
      
      pst_all=
      
      if ind == 1 
          temp=(Pa*pst')*(1../sig).*exp(-.5*( y(it)*ones(ns,1)- b'*x(it,:) ).^2../sig.^2)*sr2pfac;
          sp(it,:)=temp'; 
      elseif ind == 2
          temp=    
      end 
      
      
       %Step 3: Find the Likelihood 
       %----------------------------
       % f( y(t) | y(t-1) ) = sum( f( y(t) , s(t) | y(t-1) ) )
       
       lh(it)=sum(sp(it,:)); % i.e., lh is a weighted sum of the Gaussian likelihoods for
                            % for the separate states.
                            
       % Step 4: Update to find p( s(t) | y(t) ) 
       % ----------------------------------------
      sp(it,:)=sp(it,:)/sum(sp(it,:));
      pst=sp(it,:);
   end
end
